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Abstract 

We improve the accuracy of the extrapolation of the gluon distribution of the proton 
to very small x, and show that the charm production cross section, needed to calculate 
the cosmic ray-induced 'atmospheric' flux of ultrahigh energy prompt and Ur neutri- 
nos, may be predicted within perturbative QCD to within about a factor of three. We 
follow the sequence of interactions and decays in order to calculate the neutrino fluxes 
as a function of energy up to 10^ GeV. We also compute the prompt Uj- flux from bb 
production, hadronization and decay. New cosmic sources of neutrinos will be indicated 
if more prompt neutrinos are observed than predicted. If fewer neutrinos are observed 
than predicted, then constraints will be imposed on the nuclear composition of cosmic 
rays. The advantages of studying Vr neutrinos are emphasized. We provide a simple 
parameterization of the prediction for the inclusive cross section for c quark production 
in high energy proton-air collisions. 



^To appear in the special issue of Acta Physica Polonica to celebrate the 65th Birthday of Professor Jan 
Kwiecihski. 
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1 Introduction 



Very high energy 'cosmic' neutrinos with energies in excess of about 10 TcV offer a unique source 
of valuable information about energetic events far away in the Universe; see, for example, the 
reviews in Refs. [1, 2]. This has led to the development of neutrino telescopes, which use photo- 
multiplier tubes to detect the Cerenkov radiation emitted from the charged leptons produced 
in charged-current neutrino interactions in a large volume of water or ice, deep underground; 
see, for example, the reviews in Refs. [2, 3]. If we consider muon neutrinos, then up to about 
100 TeV the spectrum is dominated by atmospheric neutrinos from the decays of pions or kaons 
produced by cosmic ray interactions in the Earth's atmosphere. At higher energies the increased 
lifetime of these mesons means that they interact before they have the opportunity to decay. 
Above 500 TeV the decays of the much shorter-lived charmed particles become the dominant 
source of atmospheric muon neutrinos. These are known as 'prompt' neutrinos^. Their energy 
dependence follows the original cosmic ray spectrum, while the spectrum of 'conventional' 
atmospheric neutrinos falls off by an extra power of the energy because of the competition 
between the decay and interaction of the parent meson. It is clearly essential to quantify the 
flux of 'prompt' neutrinos as accurately as possible, since they provide the background to the 
sought-after 'cosmic' neutrinos. 

There exist many models of the 'prompt' neutrino flux in the relevant 10^-10^ GeV energy 
range, which yield predictions which differ by more than two orders of magnitude. Some of the 
models are purely phenomenological and have arbitrary continuation to high energy from the 
domain constrained by accelerator data. However, it has been noted that, with the inchision 
of the ncxt-to- leading order (NLO) contributions [4], perturbativc QCD gives a satisfactory 
description of the observed features of the available accelerator data on charm production, see 
for example [5]-[ll]. Moreover, the simplifled LO calculation reproduces the same behaviour 
when multiplied by an overall K factor, K ~ 2.3, to account for the NLO contribution. Several 
authors have therefore used perturbative QCD to predict the prompt neutrino flux [7]-[ll]. 

The perturbative approach, however, faces the same problem of extrapolation to high en- 
ergies. The LO diagram for forward high energy charm production is shown in Fig. 1. The 
cross section may be written in terms of the Fcynman variable xp — Pl/Pl^^i where is 
the longitudinal momentum of the charm quark; at high energies xp ~ Ec/E^ where E is the 
incident proton energy and Ec is the energy of the charm quark. Using the notation of Fig. 1, 
we have 

^{pp^ c+ ...) = / dxidx2dz g{xi,iij,)^^^^f^ g{x2,ii%)5{zxi- xf), (1) 
axp J dz 

where da /dz = sda/dt with z = (m^ — t)/s, and where g{x) is the gluon density of the proton. 
The Mandclstam variables s and i refer to the gg cc subproccss.^ The problem is that in 

^For tau neutrinos, we will see that prompt production dominates at about 10 TeV and above. 
^Throughout we take the mass of the charm quark to be rric = 1.25 GeV, following Ref. [9]. We know that 

this value, taken together with the NLO contribution {K ~ 2.3), gives a good description of all the available 
fixed-target data for cc, or rather D meson, production (which are in the region E ~ 250 GeV) [9] . 
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the high energy domain of interest we sample the gluon density at very small X2] to be specific 
X2 — MlJ^xps ~ 10~^-10~^, where \fs is the total pp cm. energy. There are no data which 
determine the gluon for x < 10^*^, and, as a rule, parton distributions for x < 10~^ are not 
available. For example, Ref. [11] shows a range of predictions for the prompt flux neutrino 
obtained by continuing the gluon distribution below x = 10^^ using the power law dependence 
xg ~ x~^ with A in the range 0-0.5. Of course, the prediction depends crucially on the value 
taken for A, and at the highest neutrino energy shown, 10^ GeV, the rates span about two 
orders of magnitude. 

P 



P 

Figure 1: The lowest-order diagram for cc production in high energy pp collisions. The 'small 
x' gluon has typical values X2 ^ M^/2sxf, where xp ^ 0.1. 

The goal of the present paper is to diminish the uncertainty in the predictions of the prompt 
neutrino flux. The major problem is to obtain the most reliable method of extrapolation based 
on the present understanding of the small x regime. In order to do this we begin, in Section 2, 
by comparing different physically-motivated extrapolations of the gluon to very small x: 

(i) DGLAP gluon with a double leading log (DLL) extrapolation, 

(ii) unified DGLAP/BFKL gluon [12] with a;^^ extrapolation, 

(iii) extrapolation with saturation effects included. 

In Section 3 we compare the predictions for the xp distribution of charm quarks produced 
in high energy pp collisions {E ~ 10^ to 10^ GeV) using the three models for extrapolating the 
gluon into the small x regime, typically x ~ 10~^ to 10~^. We argue that the extrapolation 
based on model (iii) is the most reliable, and so for the remainder of the paper we show results 
for this gluon. To determine the prompt neutrino flux we need to extend the calculation to high 
energy proton-air collisions. This is done in Section 4, where we also consider the uncertainties 
associated with the extrapolation based on model (iii). 
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In Section 5 we describe the formalism that we shall use to calculate the prompt lepton 
spectra. Starting from the production of cc pairs from the incoming cosmic ray flux, we allow 
for their fragmentation into charmed {D"^, D^, Dg and Ac) hadrons, and for their subsequent 
semileptonic decays. We include the effects of the lifetime of the charmed hadrons and, also, 
for the lifetime and decay modes of the r lepton in the Dg rUr decays. The results of the 
calculation of the prompt and fluxes are presented in Section 6, and compared with the 
'conventional' atmospheric fluxes. These latter fluxes arise from tt, K . . . decays and — u^- 
oscillations respectively. We flnd that the prompt Ur spectrum for E > 10^ GeV hes much 
above its conventional atmospheric background, whereas for the prompt i/^ spectrum this is 
only achieved for £■ > 10^ GeV. The origins of the prompt u-,- flux are the Dg — > decays 
which occur with a branching fraction of 6.4 ± 1.5% [13]. It is relevant to note that high 
energy z/^'s, unlike z/^'s, are not depleted in number by absorption in the Earth. They will 
always penetrate the Earth due to the u^- ^ t u^- . . . regeneration sequence [14]. This is 
clearly crucial for upgoing high energy neutrinos through the Earth and can be important for 
horizontal neutrinos (particularly of high energy) entering a detector deep underground. 

At first sight, it appears that the prompt neutrino fiux from bb production will, relative 
to the flux of cc origin, be suppressed, first, by a factor of order ml /ml and, second, by the 
fact that the gluon is sampled at larger x. However, for the prompt Ur fiux of bb origin, the 
suppression is partly compensated by the existence of significant TUr semileptonic decays of 
all the beauty hadrons (5=*=, Bg and A;,), in contrast to just the Dg — > tVt decays for 
charm. We calculate the prompt flux arising from bb production, fragmentation and decay in 
Section 7. Finally, in the concluding section, we discuss the implications of our results for 
neutrino astronomy and cosmic ray physics. Also, there, we summarize the uncertainties in the 
predictions of the prompt neutrino fluxes. 



2 The gluon at small x and high energy cc production 

As mentioned above, one possible method of extrapolating the gluon into the x < 10~^ regime 
is to resum the leading q;5 In In l/x terms within the DGLAP framework, which leads to a 
small X behaviour^ 



xg{x,Q'^) ~ Xog{xo,Ql) exp 



,^ln^ln^|. (2) 

\ Ois[QQ) Xq 



^In practice, we use the DLL continuation to x < 10 ^ in the form 



i I n^^ f A^iVc . as{Q) , x 16iVc , as{Q) , 10"^ 
g{x,Q ) = 5(a; = 10 ) exp W— — In — 777^ In W— — In — ^77^ In , 

where the LO coupUng as(Q) = 'iTT/{b\og{Q'^ / A^^^^)) with n/ = 4 and 6 = 25/3. Wc take Qf, = 1 GcV^ 
and xo = 0.25. We use MRST2001 partons [15] with Aqcd = 220 MeV. The rehabihty of this form of DLL 
extrapolation has been checked using GRV partons [16] which are tabulated down to a; = 10~^. 
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We denote this extrapolation by MRST on the Figures below. 

As far as we fix the scale and extrapolate to much smaller the leading contribution 
comes from q;5 In l/x terms, which, at leading order, are resummed by the BFKL equation [17]. 
Therefore a more reliable extrapolation is obtained by solving a unified DGLAP/BFKL equa- 
tion [12] for the gluon. This equation is written in terms of the gluon, unintegrated over its 
transverse momentum, which should be used with the off-mass shell^ matrix element for the 
hard gg cc subprocess amplitude [18]. In this way the result embodies the main part of the 
NLO DGLAP correction. Besides this, the unified equation includes a kinematic constraint (or 
consistency condition) [19] which accounts for the major part of the higher-order corrections to 
the BFKL evolution.^ Indeed, the power behaviour generated corresponds to A ~ 0.3 which is 
much less than the LO BFKL behaviour x~'^° with ujq — V2as\'Q.2 / . Moreover, the charm pro- 
duction cross section da/dxp calculated in terms of the KMS unintegrated gluons [12] is found 
to coincide, within 10% accuracy, with the prediction obtained with conventional DGLAP glu- 
ons [15] for 10^-10^ GeV laboratory energies, corresponding to a; = 10~^-10^'^ for which deep 
inelastic accelerator data exist. Note that, in this comparison, the prediction based on conven- 
tional partons was calculated at LO and multiplied hy K = 2.3, which, according to Ref. [9], 
accounts well for the NLO corrections. On the other hand the unified DGLAP/BFKL prescrip- 
tion already incorporates the main NLO effect at small x, and so no K factor is present in this 
X domain. We denote the results obtained using this gluon by KMS on the Figures below. 

At 10^ GeV these predictions should be reliable. However, as we proceed to higher en- 
ergies, we sample gluons in smaller and smaller x regimes with increasing gluon density, 
and so we must account for the effect of saturation. To study this effect we start with the 
Golec-Biernat, Wiisthoff (GBW) model of deep inelastic scattering (DIS) [21] (and diffractive 
DIS [22]). Let us outline the basis of the model, as apphed to qq production in DIS. The 
production of qq pairs is described by the probability of the formation of the pair by the ini- 
tial photon multiplied by the cross section for the gg-proton interaction, a. The first stage is 
given by the effective photon wave functions ^t,l (for transverse, longitudinal polarisations), 
which depend on the momentum fraction z carried by the quark and the transverse separation 
r between the q and q. The deep inelastic cross sections have the form [23] 



with X = X2; where fg is the unintegrated gluon density as defined in [12], and where (Ton, off are the gg cc 
cross sections with the small x gluon on-, off-mass-shell. The large Xi gluon, which is clearly in the DGLAP 
regime, is always taken to be on-mass-shell. 

^Explicit expressions for the next-to- leading log BFKL terms can be found in [20]. 




(3) 



^That is we use the replacement in (1) 
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For small r, the dipole cross section a is proportional to r^. To allow for multiple interactions, 
Golec-Biernat and Wiisthoff [21] parameterize a by the form 



^2 



a(x,r) = '-o^l-expj^-^JJ, (4) 



with an x-dependent saturation radius 



I f X \V2 



Qq \XqJ 

The parameterization is a simplified version of the well-known Glauber expression for, say, 
describing the multiple interactions of a pion passing through a nucleus 

<JnA = / d\ [1 - exp i-a^NTik))] . (6) 

The integral / (Pbt gives the nuclear area nR\, which is replaced by ctq in (4), and the mean 
nucleon density {T) is parameterized by l/Ro{xy, modulo normalisation. cTttjv is the 7r-nucleon 
cross section, which is equivalent to the gg-dipole cross section in the GBW model. That is, the 
exponent (Tt^nT in (6) is equivalent to {7r^r'^as/S)xg/ao, where the gluon density xg/ao plays 
the role of the mean nucleon density ( T" ) , and where the factor in brackets plays the role of 
(j^jv- It is because the gluon density grows as x decreases that we have an x dependence in the 
argument of the exponential in (4) . 

The GBW model has recently been realised [24] in terms of the gluon density, including the 
DGLAP In evolution of g{x, Q"^)- Actually in this improved form it should be considered, not 
as a model, but as a complete perturbative calculation, which in addition to the conventional 
LO collinear approach also accounts for the rescattering of the incoming qq pair. 

The power of x in (5) reflects the power growth of the gluon density in the small x region. 
The parameters ctqjXo and A were fitted to describe the small x DIS data [21]. It was shown 
that, up to ~ 20 GeV^, a good description can be achieved, even without accounting for 
DGLAP evolution. Interestingly, the value of the power, A = 0.28, turns out to be close to the 
value found in Ref. [12]. 

So far we have considered absorption for DIS. Here we are interested in gg — > cc, and not 
^g — > cc. It is therefore necessary to multiply the ^g cc cross section by 



1 - lz{l-z) 



(7) 



a el 



where the first factor is due to colour, and the second term in square brackets accounts for 
99 9 ^ cc production (see, for example, [25]). 

Note that the approach of Golec-Biernat and Wiisthoff includes only the rescattering of 
the cc pair and neglects the enhanced Reggeon diagrams which account for the rescattering of 
the more complicated Fock components of the photon (gluon) like ccg, ccgg, etc. These extra 
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components have a larger absorptive cross section. In other words, when the gluon density 
becomes sufficiently large, we must allow for gg recombination, which diminishes the rate of 
cc production. From this t channel viewpoint, the absorption is driven by the triple- Pomeron 
interaction. With the help of the Balitskii-Kovchegov equation [26, 27], we may sum up the 
resulting fan diagrams (formed from different networks of Pomeron-Pomeron recombinations 
into single Pomeron exchanges). This effect has been studied numerically in recent papers [28, 
29, 30]. However, the approach is based on LO BFKL, and does not account for the NLO 
corrections, which are known to be large. In this case we cannot simply rescale the LO prediction 
by taking a lower Pomeron intercept, cuq. The problem is that the triple- Pomeron vertex is not 
known at NLO. 

There are reasons, both from phenomenology [31] and from perturbative evolution [32], to 
believe that the triple-Pomeron coupling is small. Nevertheless, at very high energy, we expect 
the resulting absorption to be stronger. From this point of view we may regard the prediction 
based on the GBW model as the upper limit for cc production. Later, for a more realistic 
estimate of the cross section for cc production, we take account of the triple-Pomeron vertex 
by replacing (5) by^ 

Rl = ^(c+{x/xo)') (8) 

for X < 10~^ — 10"^ with c ~ 0.05-0.2. This is to protect Rq becoming too small for small x. 

The only reason why the above upper limit may be exceeded arises because the GBW 
saturation model [21] was formulated at fixed impact parameter, and so does not allow for the 
growth of the proton radius Rp with increasing energy. The radius Rp can be determined from 
the slope B of the elastic pp cross section, 

B = Bo + 2a'lnE, (9) 

where a' is the slope of the Pomeron trajectory, and E is the proton energy in the laboratory 
frame. Indeed, for a large-size dipole the GBW model saturates at cr = uo = 29 mb, whereas the 
normal soft hadronic cross sections, which should be equivalent to large-size dipoles, continue 
to grow logarithmically with energy. From a physical point of view, the normalisation ctq in (4) 
is related to the proton area ttR^. Of course, we only have the inequality (Jq < ttR^ oc B, since 
charm production originates mainly from the centre of the proton. However, since nRp grows 
with energy, a conservative upper limit is obtained by multiplying the prediction obtained from 
the GBW model by the factor B{E)/B{Eo), with Eq ~ 10^ GeV, typical of the HERA domain 
where the parameters of the model were tuned. 

A lower limit to cc production may be obtained if we assume a scaling behaviour for 
dric/dxp = {da{pp c -\- . . .) / dx) / ai^^^h where ciinci is the total inelastic pp cross section - 

^Such a form is motivated by the results of Ref. [33]. Recall that the combinatorial factor which corresponds 
to the fan diagrams is N\/N\ — > 1. Thus, contrary to the eikonal form, where we have 'Y^N{—g)^/N\ = exp{—g), 
here wc deal with a geometrical series type of expression ^j^{—g)^ = 1/(1 + ff)- Therefore we choose form (8) 
with the constant c as evaluated in [31, 32]. 
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that is if we assume that dric/dxp is independent of energy. Hence the lower hmit is 



da{E) 
dxp 



da{Eo) aineijE) 

dXp 0-inel(£'o) ' 



(10) 



normaUscd in the region E ~ 10^ GeV (x ~ 10~"^) where the parton distributions are known. 
To be more precise we should replace a-^^ci in (10) by the cross section corresponding to the 
Fock component of the proton wave function which contains charm. However, the cross section 
for this component will grow with energy faster than Uinei, and thus (10) may be regarded as an 
extreme lower limit for the charm yield. We consider the Fock charm component to be generated 
perturbatively. In principle, it would be possible to have also a non-perturbative 'intrinsic' 
charm component [34], although there is no firm experimental evidence for its existence. Such 
an intrinsic charm cross section would originate from the non-perturbative large-size domain, 
controlled by o"incij and hence its contribution would become less important, with increasing 
energy, than the perturbative cross section. 

3 Predictions for high energy cc production 

In Fig. 2 we compare the predictions for the xp distribution of charm quarks produced in pp 
collisions, as given by the three models for extrapolating the gluon^ to small x. For laboratory 
energies E ~ 10^-10^ GeV we sample the gluons in the x region 10~^-10~^ where the parton 
distributions are known from global analyses. Hence, since each model reproduces the same 
data, they give essentially the same predictions for cc production. Recall that the LO DGLAP 
result, based on MRST partons, was multiplied by a factor of 2.3. It was shown [9] that 
such a constant K factor reproduces well the NLO perturbative QCD prediction and gives 
a good description of the available fixed-target data for cc, or rather D meson, production 
for E ~ 250 GeV. Recall that, following [9], we take the mass of the charm quark to be 
rUc — 1.25 GeV. Although we use a constant K factor, K — 2.3, we have confirmed that 
the use of the parameterization of the K factor, K{Ec,Xc), given in eq. (3.4) of [8], does not 
appreciably alter any predictions. 

Up to E 10^ GeV, the GBW saturation model practically coincides with the DGLAP 
(MRST) prediction. For higher energies the GBW cross section is lower due to absorptive 
effects. On the contrary, the prediction based on KMS partons becomes higher, as well as lower 
at the lower energies. The 'unified' KMS evolution includes the BFKL Inl/x resummation and 
generates a power growth, a;"'^, of the gluon density as we extrapolate to small x. This evolution 
embodies a kinematic constraint (or consistency condition) which accounts for a major part of 
the NLO and higher-order BFKL effects. However, the power, A ~ 0.3, is appreciable, and the 
growth exceeds the double logarithmic DGLAP growth of the MRST result. Another feature 

^The gluons in (1) are evaluated at a scale /xf equal to the transverse mass of the charm quark for the MRST 
and KMS models; that is fij^ = + PcT- the GBW extrapolation we take fxp = {^/r), where r is the 
separation between the c and c quarks. 
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to note is that the shape of the xp distribution becomes a httle steeper with increasing energy, 
as seen in Fig. 2(d), which shows the predictions of the GBW extrapolation for three different 
energies. 

On the other hand, at low energies, we see from Fig. 3 that the KMS prediction falls away. 
Indeed, it is about a factor of two below the GBW/MRST predictions for Ec = 10^ GcV. 
Here, we are sampling the gluon at x values above 0.01, where the In l/x resummation is not 
effective and where LO DGLAP evolution takes over. To put it another way, it was observed 
in the KMS unified BFKL/DGLAP approach [12], that the cross section for an off-shell gluon 
is enhanced by In l/x effects for x < 0.01, whereas it rapidly tends to the LO on-shell DGLAP 
formula as x increases above this value. Hence in the low energy regime a factor X ~ 2 should 
be included in the KMS prediction, with the factor dying away with increasing energy, as we 
enter the BFKL regime. 

A convenient way to summarize the relevant energy behaviour of the da'^/ dxp cross section 
is to plot the moment' [35] of the xp distribution, see for example [36]. For high energies 
{E > 10^ GeV) the incoming primary cosmic ray flux falls down as with 7 = 2.02. 

Therefore the charm flux is proportional to the moment 



This moment is shown in Fig. 3, where the difference between the saturation model and the 
other two models becomes apparent for Ec > 10^ GeV. Note that here we flx the energy Ec of 
the outgoing c quark, rather than that of the incident proton which was used in Fig. 2. 

Although the GBW model predicts the smallest cross section of the three models, it should 
be considered as the upper limif for cc production as it only accounts for part of the absorptive 
effects. For the reasons mentioned earlier, the GBW model is more than a model - rather it 
should be regarded as a ftill leading-order QCD prediction with cc rescattering effects included. 
Of course there is, in addition, absorption of the gluons in the evolution process. It appears 
likely that the consequent reduction of the cc cross section due to this additional absorption 
may be partially compensated by the growth of the proton radius with increasing energy. We 
investigate this, and other effects, in the next Section; see Fig. 4 later. Therefore, from now 
on, we will base our study on the GBW approach and its variations. 

4 cc production in proton— air collisions 

For a precise study we need the charm yield in p-nuclear (air) collisions. An advantage of 
the GBW saturation model is that it may be straightforwardly extended from pp to p-nuclear 

^Modulo a possible growth of the proton radius, see (9). 




(11) 
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collisions. Within the GBW framework we may account for the eikonal rescattering of the cc 
pair within the nucleus by replacing d in (4) by 



oo 



(-(A-l)(T/87rS^) 



n-l 



nn\ 



(12) 



n=l 



where, for air, the mean atomic number A = 14.5 and the slope Ba = (r^)/6 = 29 GeV~^. 
Note that in the numerator we have {A — 1) and not A. In this way we exclude rescattering 
on the nucleon where the cc pair is created. This rescattering is controlled by a different slope 
(t^ Ba), and is already accounted for in the cross section given in (4). We have taken the root- 



mean-square nuclear radius^ y (r^) = 2.6 fm, and assumed a gaussian distribution of the nucleon 
density for light nuclei. Note that the replacement of a occurs before the integration over the c-c 
separation in (3). In summary, the inclusive cc cross section for proton-air interactions is given 
by the sum of proton-nucleon cross sections, with the only nuclear effect being the enhanced 
absorption of the produced cc pair. In the Appendix we give a simple parameterization which 
reproduces the proton-air cc production cross section in the energy interval 10^ < E < 10^^ GeV 
to within ±5%. 

The 7 = 2.02 moment of charm production in high energy proton-air collisions is shown 
in Fig. 4. The dashed curve is the prediction of the extrapolation based on the original GBW 
saturation model, while the upper dotted curve (marked by B) includes the possible effect of the 
growth of the proton radius with energy, as discussed at the end of Section 2, see in particular 
(9). Instead of using the Pomeron slope, a' — 0.25 GeV~^, measured in elastic pp scattering, 
in (9) we use the value 0.11 GeV~^ together with slope B = 4.4 GeV~^ at = 90 GeV, which 
were deduced from elastic J/ t/j production data at HERA [38] . These values are appropriate for 
charm production. The lower dotted hue (marked gsp) reflects the inclusion of an additional 
absorptive effect — the absorption of gluons. We take c = 0.2 in (8), which corresponds to the 
largest estimate^°, g^p = 2 GeV~^, of the triple-Pomeron vertex. This choice is made to show 
the full extent of the uncertainty in cc production. If both the above effects are included (the 
radius growth and g^p), then the solid curve (marked B + g^p) is obtained. Conservatively, we 
predict that will lie within the shaded region in Fig. 4; the most likely behaviour is that it 
will follow the cross-hatched region. Even the conservative prediction has much less uncertainty 
than previous estimates. However, for completeness, the dot-dashed curve shows the extreme 
lower limit, given by scaling formula (10), but where now Uinei is the proton-air cross section. 
All the variations of the original GBW model were normalized in the region E ~ 10^ GeV, 
where the partons sampled in the hard subprocess are known, and the model tuned to the data. 

^For example, Ref. [37] gives the r.m.s. radius of oxygen as 2.7 fm. To obtain the result for air we take the 
usual r oc ^3 dependence. 

^"Triple- Pomeron phenomenology gives g2.p ~ 0.5 GeV~^. However, the hare triple-Pomeron vertex may be 
a few times larger, since the phenomenological value already accounts for some screening effects. It was argued 
in Ref. [32] that in perturbative QCD we expect gzp = (0.5-2) GeV~^. 
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5 Prompt neutrinos: development of the air shower 



Our aim is to predict the spectra of prompt and neutrinos produced in the atmosphere by 
cosmic rays. Prompt leptons originate from the following sequence: the production of cc pairs, 
their fragmentation into charm hadrons which then undergo semi-leptonic decay. In the lower 
energy range {E < 10'' GeV) it is possible to estimate the Icptonic spectra by simply taking 
a product of moments of the various distributions, see, for example, Ref. [36, 39]. However, 
for E > 10^ GeV the decay length of D mesons becomes comparable with the depth of the 
atmosphere, and so it is necessary to follow the development of the air shower in more detail. 
It is described by a set of equations in terms of the 'depth' X of the atmosphere transversed 
by a particle in the vertical direction. X is related to the height h by 



X 



POO 

/ p{h')dh' (13) 

J h 



where p{h) is the density of the atmosphere at vertical height h. We take the same exponential 
profile of the atmosphere^ ^ as was used in Ref. [7]. The sequence of equations determine 
(paiE, X), which are the fluxes of the corresponding particles with energy E at depth X, where 
a = N,c,i,l (that is nucleon, c quark, charmed hadron, lepton). The initial flux (f)i^{E,0) is 
the known primary cosmic ray flux. All other initial fluxes are zero, that is (paiE, 0) = for 
a = c,i,l. The set of equations which determine 4>n —>■ 4>c 4>i 4>i are 

- ' :ME,xy. (14) 



dX An{E) 



,{E,X) = / dE'dx^ME'.X)'^^ 5{E-x,E')- (15) 
JE p aXc 



dX Ai{E, X) 

with i = D^, L>°, Df, A^; and finally 



POO rlrj ^ 

{E,X)+ dE'dx(l)c{E',X)^ S{E-xE') (16) 
Je dx 



{E,X) 



dX 



^J-4E'a.ME',X)^^^^B(i D'^'^E - .E') (17) 



where B{i ^ I) is the branching fraction of the decay of the charmed hadron i to lepton I. The 
nucleon attenuation length is 

A.W-^^, (18) 

J- — ^NN 

where Z^n is the spectrum- weighted moment for nucleon regeneration and Aat is the interaction 
thickness 

Aiv = p{h) / Y.^NA{E)nA{h). (19) 

/ A 



^^Wc have checked that the numerical results essentially do not depend on the precise parameterization of 
the density of the atmosphere. 
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nA{h) is the number density of air nuclei of atomic number A at height h and a^A is total NA 

inelastic cross section. Instead of the sum over A in (19), we take the mean value (^4) = 14.5 
for air. Note that the factor ua/p in (15) arises from (19). For the nucleon-air cross section, 
(TAT-air, wc take the parameterization of Bugaev et al. [40]. For the incoming cosmic ray flux 
we take the parametrisation given in [36] denoted as TIG with knee. Also from [36] we take 
parametrisation of Zj^^ which depends on energy and takes into account the knee, which is 
consistent with [7]. This Z factor includes the regeneration by the p,n,N* . . . particles. 

The attenuation length Ai{E,X) of the charmed hadrons consists of two parts: the decay 
length Xf^'^ and attenuation due to their strong interactions with air nuclei. The decay lengths, 
Xf^, of the various charmed hadrons depend, via their X dependence, on the density of the 
atmosphere, 

Af- = cr,^p(X), (20) 

where rrii and Tj, the mass and lifetime of the i^^ charmed hadron, are taken from [13]. The 
attenuation due to the strong interactions of the produced charmed hadrons with the air has a 
form similar to (18), namely Ai/(1 — Zee)- To calculate we assume an absorptive cross section 
equal to half the absorptive p-air cross section (based on quark counting) , and we take a charm 
regeneration factor Z^c — (O.S)'''. That is, we estimate that the leading charm quark carries a 
fraction x — md {rric + mg) ~ 0.8 of the incoming energy, where ruq ~ 0.3 GeV is the mass of 
a light constituent quark. For simplicity, we take the same Zee and Aj for all charm hadrons; 
for Ac we expect the larger cross section to be approximately compensated by the larger Zee- 
Thus, finally, Aj is given by 

I I I - Zee 

From (14) the light baryon flux is given by 

ME.X) = 0Ar(£;,O)exp(-X/Ajv(£;)). (22) 
If we insert (15) into (16), then the individual charm hadron fluxes are 

I dX" \ 

ME, X) = I dX' exp (- ^^^j S,^,(E, X') (23) 

where 



S^^,{E,X)= dE'^^^ '^'^^ME',X) (24) 
Je p E' ax 

and X = E/E'. 

The spectra of charmed hadrons da^^'^/dx are calculated using the three models shown in 
Fig. 2. The ratios of the different charm yields (after the hadronization of the c quark) are 
given in Ref. [6] to be 

a(D+,D0)=°-2' ^^^t;^ = 0.3, ^-0.5. (25) 
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More recent data [41] favour a smaller value^^ of the first ratio quoted in (25), namely 

(-) 

We take this value in our analysis. 

Note that the baryon is produced by the recombination of a c quark with a spectator 
diquark of the incoming nucleon^^. It is not produced from a c quark. Therefore the parton 
momentum fraction x\ carried by the Ac is x\ = X(i + x^. The diquark momentum fraction x^ 
will be less than |(1 — Xc), as part of the energy is carried away by the third valence quark (the 
factor |), by the c quark and by gluons. To allow for this, we therefore take = ^(1 — Xc), 
which leads to 

XA = |(1 + Xe) (27) 

This is found to be in good agreement with the distribution generated by the PYTHIA Monte 
Carlo [42], which has a maximum in the region xa — 0.5-0.6. Of course a very slow c quark 
is unlikely to combine with a fast light diquark (with Xd ~ 0.5). Therefore we introduce a 
cut-off, Xc > in (27). Assuming a mean velocity of the c quark to be (t)^) ~ 0.25, we 
estimate x° ~ 0.1. For A;, production, which we discuss below, the heavier b quark will carry 
a larger fraction of the A;, momentum. In this case we take the cut-off to be = 0.25. 
As a check, we also compute the prompt flux, arising from c — > A^, using an alternative to 
prescription (27). We assume that for the diquark x^ = {md/rric)/ Xc with a constituent diquark 
of mass rud = 2mq = 0.7 GeV and rric = 1.5 GeV. With this assumption 

XA = lA7xc, (28) 

or, in the case of beauty, xa = I.I6X5, using nib = 4.5 GeV. 

Rather than using a fragmentation function for drf^^ / dx, for D mesons we take Xd = 
0.75a;c. This is sufficient for our purposes^^. For illustration we compare, in Fig. 5, the prompt 

+ flux at ground level, 0^^_|_p^(i?), obtained from (17) using different forms of fragmen- 
tation. Clearly the upper curve, corresponding to no fragmentation, gives an overestimate of 
the flux. Moreover, due to the presence of additional hght "sea" quarks, we expect a harder 
distribution for the fragmentation in pp production than that obtained in e^e~ collisions (lower 
curve). Hence our use of fragmentation corresponding to one of the middle curves. We see that 
both the Ac hadronization prescriptions (27) and (28) give very similar fluxes. At the highest 
energies the fraction of neutrinos coming from Ac (relative to those from D) increases due to 

-'^^For isolated c-quark fragmentation, PYTHIA [42] gives a ratio of 0.08. 

^^The PYTHIA/JETSET Monte Carlo gives only 2.5% Ac baryons in the fragmentation of an isolated c jet, 
that is in the absence of spectator diquarks. 

^^That is the distributions dn''^'^/dx were taken to be proportional to 6{xd — O.TSXc) for D mesons, whereas 

for Ac we assume that they are proportional to 5{xa — ^{xc + 1)) for Xc > 0.1, or S{xa — 1.47.Tc). We note 
that PYTHIA gives a harder x = Ed / Ec distribution than that shown for the Peterson et al. function [43] in 
Ref. [13]. 
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the short Ac hfetime. Therefore the choice xa — ^{1 + Xc) of (27), which corresponds to a larger 
(xa) than (28), gives a larger neutrino flux for Ei, > 10^ GeV. The results below correspond 
to using prescription (27). 

For the leptonic decay of each charm hadron i, the distribution dn^'^^ /dxi was generated 
by the PYTHIA Monte Carlo [42]. Note that in PYTHIA 6.2 the Ds TVr branching ratio 
was set to be 1%, whereas the latest value is 6.4 ± 1.5% [13]. Since Ds tv^ is almost the 
only source of prompt neutrinos, it is important to renormalise the yield using the updated 
branching ratio. It is interesting to note that the -D+ decay produces more prompt neutrinos 
than z/^, since the from r"*" decay has the large xi. Of course, the reverse is true for D~ 
decay. 

6 Prompt Vy, and fluxes 

We present the predicted yields of prompt z/^ and Vr neutrinos from charmed hadrons produced 
by cosmic ray collisions in the atmosphere. The flux of z/g neutrinos is essentially equal to that 
of z/^ neutrinos. High energy prompt electrons are completely degraded in the electromagnetic 
cascade. For prompt muons the electromagnetic interaction is much weaker; it was demon- 
strated in Ref. [44] that the prompt \x flux is only about 10% smaller than the prompt flux 
at the surface of the earth. 

In Fig. 6 we plot the prompt z/^ and fluxes predicted by the three models for the ex- 
trapolation of the gluon distribution to very small x. Although in Figs. 3 and 4 we compared 
different models of extrapolation using a flxed 7 = 2.02, in the actual calculation of the neutrino 
fluxes we used the observed primary cosmic ray flux, which corresponds to different values of 7 
above and below the 'knee'. The sharp fall-off for neutrino energies E > 10^ GeV is due to the 
increase in the decay length of the charmed hadrons, arising from Lorentz time dilation. Clearly 
in the gluon kinematic domain of interest (x < 10^^ and < 10 GeV^) saturation efi^ects be- 
come important. The reason for the behaviour of the KMS prediction — below at low energies, 
above at high energies — was given in the discussion concerning Fig. 3. We have argued that 
extrapolations based on the GBW model and its variations, as shown in Fig. 4, give the most 
reliable predictions for £■ > 10^ GeV. Thus in Fig. 7 we show the spread of predictions of the 
neutrino fluxes based on the shaded domain in Fig. 4. In Fig. 7 we also show the conventional 
atmospheric flux of (from tt and K decays, etc.). Moreover, there is a small probability that 
atmospheric z/^ neutrinos may oscillate into Ur neutrinos and so provide an 'atmospheric' z/,- 
flux. We calculate this flux using the 3a ranges of the sin^ 6'atm and Am^rpj^ neutrino mixing 
parameters found in an analysis [45] of the Super-Kamiokande [46] and MACRO [47] data. The 
resulting atmospheric z/^ flux is shown by a band in the lower plot of Fig. 7. 

We discuss the neutrino fluxes of Fig. 7 in the concluding section, after we have included the 
contribution to the prompt u^- spectrum arising from bb production, fragmentation and decay. 
However, first, we show in Fig. 8 the effect of charmed hadron interactions with the atmosphere. 
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The effect is illustrated by the difference between the dashed curve, for which the interactions 
are suppressed (that is Aj = 0), and the continuous curve with the interactions present. For 
this comparison we use the GBW extrapolation, the one given by the solid curve in Fig. 4. 

We emphasize that the predictions for the prompt neutrino fluxes depend strongly on the 
nuclear content of the primary cosmic rays. So far we have assumed that the cosmic rays are 
composed entirely of protons. Suppose that the protons were to be replaced by nuclei of the 
same energy E and of atomic number A. Then we have to scale the energy of the primary 
interacting nucleon to E/A. Roughly speaking, this reduces the neutrino flux 0j, in the plateau 
region (10^-10^ GeV) of the E^cpi, plot by A^. On the other hand, the number of incoming 
nucleons is A, so as a consequence we expect an A^ suppression. A detailed calculation for 
incoming nuclei with (A) — 7 [48] gives the dotted curve in Fig. 8. The suppression of the 
original GBW continuous curve (for A — 1) is apparent. 

7 Prompt flux from bb production 

At first sight it appears that if we also allow for bb production, with the same cross section, 
fragmentation and decay then we will approximately double the i^t flux. However note that, 
in comparison with charm production, the cross section for beauty production is about 30 
times smaller due to the factor ml /ml and to the larger value of X2 of the gluon structure 
function which is sampled, see eq. (1), see Fig. 9. Nevertheless, the high energy production of 
Vt neutrinos from beauty decays is not negligible. For charm, only the decay of the Dg may 
produce z/^ neutrinos, whereas now B^, B^, Bg and A;, semileptonic decays also give rise to a 
significant z/^ fiux. Indeed including B and A5 decays enlarges the predicted total prompt z/^ 
fiux by about 40% for E,, ~ 10^ GeV, and even more at higher energies, see Fig. 10(b). This is 
considerably larger than the fiux calculated in [49], where the beauty induced contribution 
grows from 1% at E^, = 10^ GeV to about 10% at Ej, = 10® GeV. The reasons why we obtain 
a larger fraction of neutrinos from beauty are, first, that we use (26), rather than (25), for 
c ^ Dg hadronization and, second, our cross section for charm production is more suppressed 
at high energies by absorptive corrections, than the more compact 66 production process. 

Fig. 10(a) shows the breakdown of the total prompt Vt flux. We take the ratios of the 
B^, Bg and A;, beauty yields (after hadronization of the h quark) to be given by exactly 
analogous relations to (25) and (26) for charm. Recall that i/^-'s of cc and hb origin come, 
respectively, from Dg tv^ and from the semileptonic tv^- decay modes of S^, B'^ ^ Bg and 
Af,. Thus we have a direct component and an indirect component coming from the r ^ v-y 
decays. For charm the former component is very small since the direct carries away a small 
energy fraction. Thus the t ^ decay gives the dominant component, which is truncated 
at an energy when the r has insufficient time to decay. For beauty the direct and indirect 
components are comparable until the energy regime is entered at which the r has not enough 
time to decay. 
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To calculate the b quark cross section we have used the same GBW model as described in 
Sections 2-4 with a quark mass rrij, = 4.5 GeV. We have approximated the b ^ B fragmentation 
function by a delta function at a; = Eb/E}, — 0.9. Also we have assumed that the hadronization 
and the cross-section of the beauty hadron-air interaction do not depend appreciably on the 
nature of the heavy quark. So we have used the same attenuation lengths Aj and the same 
ratios of the different beauty hadron yields (eqs. (25) and (26)), as were used for the charm 
case. 

We neglect the bb contribution to the (and i^e) fluxes. They never exceed about 2-3% of 
that arising from cc production and decay. 

8 Concluding discussion 
8.1 Implications for Physics 

The flux shown in Fig. 7(a), and the z/^ flux shown in Fig. 10(b), have important implications 
for neutrino astronomy. These neutrino fluxes arise from cosmic ray interactions with the 
atmosphere. They therefore provide the background to searches for cosmic neutrinos (for which 
there are many exciting New Physics scenarios). They also provide a possibility to calibrate the 
detectors of the new neutrino telescopes. A particularly interesting energy domain is where the 
prompt neutrino flux has emerged from the sharply falling 'conventional' atmospheric neutrino 
flux. Here 1/^- appears to have an enormous advantage in searching for a signal for New Physics, 
see Figs. 7 and 10. Due to mixing, the fluxes of and are equal for incoming cosmic 
neutrinos. The flux of prompt z/^ neutrinos (which arises mainly from Dg rVr decays, but 
with a significant 40% contribution from 5, fi^, Af, semileptonic decays) is about ten times 
less than that for prompt z/^, and the number of neutrinos produced from the conventional 
atmospheric flux (via —>■ Ur oscillations) is negligible for E > 10"^ GeV. There is only a 
little time for neutrino mixing in the atmosphere. Thus the tau neutrino flux originating from 
atmospheric neutrinos is greatly suppressed. As a consequence, tau neutrinos offer an ideal 
means of identifying neutrinos of cosmic origin, and for searching for New Physics. 

There is another reason to concentrate on the Ut flux. For values of E^, above 10^ GeV the 
flux is significantly depleted in the passage of the neutrinos through the Earth. The absorption 
increases rapidly as Ei, is increased [50]. On the other hand, high energy i/j- neutrinos have the 
unique advantage that they are not depleted in number no matter how much of the Earth that 
they pass through. If they suffer a charged-current interaction, they produce a r lepton which 
subsequently decays, regenerating a u^- neutrino with degraded energy [14]. It is, moreover, 
interesting to note that tau neutrinos of energies about 10^ GeV may produce spectacular 
distinctive signatures in the planned 1 km^ high energy neutrino detectors made of strings of 
photomultipher tubes distributed throughout a naturally occurring Cerenkov medium, such as 
water or ice, deep in the ocean or ice cap. A charged-current i/t interaction may produce a 
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contained 'double bang' signature [51] or a 'lollipop' signature [52]. In the former, the first 
bang corresponds to the hadronic shower produced along with the subsequent r lepton and the 
second bang is the shower associated with the r decay. A lollipop event is when only the second 
shower occurs within the detector and the r lepton which initiates the shower is identified by 
the relatively weak ionization that it causes. 

Another interesting possibility follows from the observation that electron neutrinos may 
be distinguished experimentally from muon and tau neutrinos. Electron neutrinos give rise 
to distinctive showers, which result from their charged- current interactions, which character- 
istically are contained within the large detector volume, and hence identified. The conven- 
tional atmospheric fiux of Ue neutrinos comes mostly from kaon decays and so the relative fiux 
^e/^/i ^ 0.1 [53]. On the other hand, the prompt and fiuxes are practically equal to each 
other. Hence the prompt fiux should be more visible over its conventional background, than 
the prompt Vr flux^^. Of course, the background reduction is much less than for Vr neutrinos, 
but this will be partially offset since it is likely that neutrinos will be easier to identify than 
Vj. neutrinos. 

So far we have shown the prompt neutrino fluxes arising from vertically-incident cosmic 
rays. However, the predicted fluxes are sensitive to the zenith angle. In particular, as we 
depart from the vertical incident direction, more atmosphere is encountered, which gives more 
time for the charmed hadrons to decay and hence allows more high energy neutrinos to be 
produced. This is well illustrated by Fig. 11, which shows the prompt neutrino fiuxes produced 
by horizontally-incident cosmic rays. Indeed, for E^, > 10^ GeV, the fiux in the horizontal 
direction is noticeably larger. The observation of the zenith angle dependence will allow the 
'atmospheric' background neutrinos to be identified and hence cosmic neutrinos (and therefore 
New Physics) to be isolated. 

8.2 Uncertainties in the predicted prompt and fluxes 

We have emphasized the importance to neutrino astronomy of reliable predictions of the prompt 
muon and tau neutrino fluxes, which arise from cc and bb production, hadronization and decay 
in the atmosphere. They provide the background to the search for cosmic neutrinos originating 
from Active Galactic Nuclei or elsewhere. We have argued that the predictions based on 
perturbative QCD, given in Figs. 7(a) and 10(b), have much less uncertainty than those that 
already exist in the literature. How reliable is perturbative QCD and the extrapolation of 
the gluon into the very small x domain? Perturbative QCD is expected to be valid for bb 
production^^, and should also hold for cc production since the factorization scale //f ~ ttIc- 

^^Wc thank John Bcacom for drawing this possibihty to our attention. 

^^Unlike cc production, there are no accelerator bb data in the fragmentation region which is relevant to this 
study. However, there exist Tevatron data for central bb production with, typically, pbr > 5 GeV [54] which 
appear to exceed the NLO QCD prediction. Recently it has been shown [55] that the discrepancy may be 
reduced to an acceptable level by using a fragmentation function with (x) larger than that of the conventional 
Peterson et al. parameterization [43] . 
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It is important to note that at fixed target energies, E ~ 250 GeV, the NLO predictions for 

da/dxp are in agreement with the available data, while for very high energies the mean scale, 
(1/r), increases due to saturation effects. Since the gluon distribution is determined by HERA 
deep inelastic data down to x ~ 10^"^, predictions, which use partonic structure determined 
from these data and which agree with fixed target charm data, should be reliable up to about 
E ~ 10^ GeV, and perhaps an order of magnitude or so above. Fig. 3 shows examples of three 
very different models which illustrate that this is indeed the case. (Recall that in Section 3 
we explained why the KMS prediction is not applicable at the lower energies shown in Fig. 3.) 
However, as we proceed to higher energies, and sample smaller and smaller x values and hence 
increasing gluon density, we must include the effects of saturation. Though we have labelled 
the predictions by GBW, we have not simply used the saturation model of Golec-Biernat and 
Wiisthoff, but rather we investigated the uncertainties in including saturation effects in some 
detail, see Fig. 4. We concluded that the most reliable predictions would be obtained using 
the continuous curve in Fig. 4, with conservative errors shown by the shaded band, which are 
reflected in the shaded bands on the prompt neutrino fluxes shown in Fig. 7. These bands 
represent the uncertainty in cc production. 

In addition, there are also uncertainties associated with the fragmentation of the charm 
quark (see Fig. 5) and with the D meson attenuation due to its strong interaction with the 
atmosphere (see Fig. 8). Another source of uncertainty is the x distribution of the Ac hyperons. 
The recombination of the c quark with the ud diquark of the incoming proton gives the Ac 
a rather large x. We use the simplified formula x\ — |(1 + Xc), where a fast diquark with 
Xud ~ 0.5 recombines with any c quark with Xc > 0.1. The relatively large uncertainty arises 
because the flux is proportional to the x"'^'^'^'^ moment, (11), of the x distribution. Thus 
the larger value of x carried by the Ac (as compared to the D mesons) compensates for the 
small probability of Ac formation, see (25). As a consequence, the contribution of the Ac to the 
prompt z/^ flux increases from about 10% to 40% as the neutrino energy increase from 10^ to 
10^ GeV. Unfortunately the x distribution of Ac hyperons at high energies {E > 10 TeV) is not 
measured yet. However, the x^-^^ moment given by PYTHIA is in agreement with the above 
prescription, (27), that x\ — |(1 + Xc) for Xc > 0.1. 

Overall, we see that the uncertainty in the prompt neutrino flux predictions is about a 
factor of 3. For the Vr flux from charm, which relies entirely on the Ds — > rur decay, there is 
also an extra uncertainty associated with the branching ratio and with the production ratio of 
(26). There is less uncertainty in the flux from beauty, since it originates more uniformly 
from B^, B^, Bg and A^ semileptonic TVr decays. 

Finally we note that the predictions have been made assuming that the incoming cosmic 
rays are predominantly composed of protons. If, however, it should happen that the observed 
and Vt- prompt fluxes lie below the predictions of Figs. 7(a) and 10(b), then the measurements 
may give important constraints on the nuclear composition of cosmic rays. 
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Appendix: charm production in proton— air collisions 

Here we present a simple parameterization which may be used to reproduce the inclusive cross 
section of c-quark production in proton-air collisions, based on the GBW model (which yields 
the continuous curves in Figs. 2-5). To ±5% accuracy for 0.05 < x < 0.6 we have 

x—^ ' ^ Axi^ il - x^-^Y (29) 

ax 

where [3 = 0.05 - 0.016 ln(E/10 TeV) and 

n = 7.6 + 0.025 ln(E/10 TeV) 1 „ 



A = 140+ (llln(£;/0.1 TeV))i-65 /xb 



for 10* < < 10" GeV, 



n = 7.6 + 0.012 ln(i?/10 TeV) 1 f^, ^ F < m^i TeV 

A = 4100 + 2451n(E/10^TeV)^b r°'^° 

For E = 10^^ GeV, the parameterization overestimates the cross section by about 10% for 
X < 0.2. 



References 

[1] R.J. Protheroe, Nucl. Phys. Proc. Suppl. 77 (1999) 456; 
R. Ghandi, Nucl. Phys. Proc. Suppl. 91 (2000) 453; 

J.G. Learned and K. Mannheim, Ann. Rev. Nucl. Part. Sci. 50 (2000) 679; 
F. Halzen, astro-ph/0301143. 

[2] F. Halzen and D. Hooper, Rept. Prog. Phys. 65 (2002) 1025. 

[3] C. Spicring, Prog. Part. Nucl. Phys. 48 (2002) 43. 

[4] P. Nason, S. Dawson and R.K. EUis, Nucl. Phys. B327 (1989) 49, erratum, ibid. B335 
(1990) 260. 



19 



[5] J.A. Appel, Ann. Rev. Nucl. Part. Sci. 42 (1992) 367. 

[6] S. Frixione, M.L. Mangano, P. Nason and G. Ridolfi, Nucl. Phys. B431 (1994) 453; 

Adv. Ser. Direct. High Energy Pliys. 15 (1998) 609 in Heavy Flavours II ed. by A.J. Buras 
and M. LmcJner (World Scientific) [hep-ph/9702287]. 

[7] M. Thunman, G. Ingelman and P. Gondolo, Astropart. Phys. 5 (1996) 309. 

[8] L. Pasquah, M.H. Reno and I. Sarcevic, Phys. Rev. D59 (1999) 034020. 

[9] G. Gelmini, P. Gondolo and G. Varieschi, Phys. Rev. D61 (2000) 036005. 

[10] G. Gelmini, P. Gondolo and G. Varieschi, Phys. Rev. D61 (2000) 056011. 

[11] G. Gelmini, P. Gondolo and G. Varieschi, Phys. Rev. D63 (2001) 036006. 

[12] J. Kwiccihski, A.D. Martin and A.M. Stasto, Phys. Rev. D56 (1997) 3991; 

J. Kwiecihski, A.D. Martin and A.M. Stasto, Acta Phys. Polon. B28 (1997) 2577. 

[13] Review of Particle Properties: K. Hagiwara et al., Phys. Rev. D66 (2002) 1. 

[14] F. Halzen and D. Sahzbcrg, Phys. Rev. Lett. 81 (1998) 4305; 
T. Stanev, Phys. Rev. Lett. 83 (1999) 5427; 

S.I. Dutta, M.H. Reno and 1. Sarcevic, Phys. Rev. D62 (2000) 123001; 
J.F. Beacom, P. Grotty and E.W. Kolb, Phys. Rev. D66 (2002) 021302. 

[15] A.D. Martin, R.G. Roberts, W. J. Stirling and R.S. Thorne, Eur. Phys. J. C23 (2002) 73. 

[16] M. Gliick, E. Reya and A. Vogt, Eur. Phys. J. C5 (1998) 468. 

[17] L. N. Lipatov, Sov. J. Nucl. Phys. 23 (1976) 338 ; 

E. A. Kuraev, L. N. Lipatov, V. S. Fadin, Sov. Phys. JETP 45 (1977) 199; 
I. I. Balitsky, L. N. Lipatov, Sov. J. Nucl. Phys. 28 (1978) 338. 

[18] S. Catani, M. Ciafaloni and F. Hautmann, Nucl. Phys. B366 (1991) 135. 

[19] M. Ciafaloni, Nucl. Phys. B296 (1988) 49; 

J. Kwiecihski, A. D. Martin, P. J. Sutton, Z. Phys. C71 (1996) 585; 

B. Andersson, G. Gustafson, J. Samuelsson, Nucl. Phys. B467 (1996) 443. 

[20] V. S. Fadin, M. I. Kotsky, R. Fiore, Phys. Lett. B359 (1995) 181; 
V. S. Fadin, M. I. Kotsky, L. N. Lipatov, hep-ph/9704267; 
V. S. Fadin, R. Fiore, A. Flachi, M. I. Kotsky, Phys. Lett. B422 (1998) 287; 
V. S. Fadin, L. N. Lipatov, Phys. Lett. B429 (1998) 127; 

G. Camici, M. Ciafaloni, Phys. Lett. B386 (1996) 341; Phys. Lett. B412 (1997) 396 
[Erratum, ibid. B417 (1997) 390]; Phys. Lett. B430 (1998) 349. 



20 



[21] K.J. Golec-Biernat and M. Wiisthoff, Phys. Rev. D59 (1999) 014017. 

[22] K.J. Golec-Biernat and M. Wiisthoff, Phys. Rev. D60 (1999) 114023. 

[23] N.N. Nikolaev and E.G. Zakharov, Z. Phys. C49 (1991) 607; ibid. C53 (1992) 331. 

[24] J. Bartels, K.J. Golec-Biernat and H. Kowalski, Phys. Rev. D66 (2002) 014001. 

[25] B.L. Combridge and C.J. Maxwell, Nucl. Phys. B239 (1984) 429. 

[26] I.I. Balitskii, Nucl. Phys. B463 (1996) 99; Phys. Lett. B518 (2001) 235. 

[27] Yu. V. Kovchegov, Phys. Rev. D60 (1999) 034008. 

[28] K.J. Golec-Biernat, L. Motyka and A.M. Stasto, Phys. Rev. D65 (2000) 074037. 

[29] M. Braun, Eur. Phys. J. C16 (2000) 337. 

N. Armesto and M.A. Braun, Eur. Phys. J. C20 (2001) 517. 

[30] E. Levin and M. Lubhnsky, Nucl. Phys. A696 (2001) 833; Eur. Phys. J. C22 (2002) 647; 
M. Lublinsky, Eur. Phys. J. C21 (2001) 513. 

[31] A.B. Kaidalov, Phys. Rep. 50 (1979) 157. 

[32] J. Bartels, M.G. Ryskin and G.P. Vacca, Eur. Phys. J. C27 (2003) 101. 

[33] Ya.I. Azimov, V.A. Khoze, E.M. Levin and M.G. Ryskin, Nucl. Phys. B89 (1975) 508. 

[34] S.J. Brodsky, P. Hoyer, C. Peterson and N. Sakai, Phys. Lett. B93 (1980) 451-455. 

[35] T.K. Gaisser, Cosmic rays and particle physics, Cambridge Univ. Press, 1992. 

[36] C.G.S. Costa, Astropart. Phys. 16 (2001) 193. 

[37] M. Anguiano and G. Co', J. Phys. G27 (2001) 2109. 

[38] ZEUS Collaboration: S. Chekanov et al., Eur. Phys. J. C24 (2002) 345; 
HI Collaboration: C. Adloff et al., Phys. Lett. B483 (2000) 23. 

[39] C.G.S. Costa, F. Halzen and C. Salles, Phys. Rev. D66 (2002) 113002. 

[40] E.V. Bugaev et al, Phys. Rev. D58 (1998) 054001. 

[41] E769 Collaboration: G. Alves et al., Phys. Rev. Lett. 77 (1996) 2392; 
WA92 Collaboration: M. Adamovich et al., Nucl. Phys. B495 (1997) 3. 

[42] T. Sjostrand et al., Comput. Phys. Commun. 135 (2001) 238; 
T. Sjostrand, Comput. Phys. Commun. 82 (1994) 74. 



21 



[43] C. Peterson, D. Schlatter, I. Schmitt and Peter M. Zerwas, Phys. Rev. D27 (1983) 105. 

[44] G. Gelmini, P. Gondolo and G. Varieschi, Phys. Rev. D67 (2003) 017301. 

[45] M. Maltoni, T. Schwertz, M.A. Tortola and J.F.W. Valle, Phys. Rev. D67 (2003) 013011. 

[46] Super-Kamiokande Collaboration: Y. Pukuda et al., Phys. Rev. Lett. 81 (1998) 1562. 

[47] MACRO CoUaboration: A. Surdo, talk at TAUP 2001, Gran Sasso, Italy, 
http: //www.lngs. infn.it/. 

[48] C.G.S. Costa, F. Halzen, J. Bellandi and C. Salles, Phys. Rev. D54 (1996) 5558. 

[49] L. Pasquah, M.H. Reno, Phys. Rev. D59 (1999) 093003. 

[50] R. Gandhi, C. Quigg, M. Reno and I. Sarcevic, Phys. Rev. D58 (1998) 093009; 
J. Kwiecihski, A.D. Martin and A.M. Stasto, Phys. Rev. D59 (1999) 093002. 

[51] J.G. Learned and S. Pakvasa, Astropart. Phys 3 (1995) 267; 

H. Athar, G.Parente and E. Zas, Phys. Rev. D62 (2000) 093010. 

[52] F. Halzen and D. Hooper, in preparation; see also the Ice Cube design document, 
www.icecube.wisc.edu. 

[53] T.K. Gaisser and M. Honda, Ann. Rev. Nucl. Part. Sci. 52 (2002) 153. 

[54] CDF Collaboration: F. Abe et al., Phys. Rev. Lett. 71 (1993) 500; ibid. 71 (1993) 2396; 
ibid. 75 (1995) 1451; D. Acosta et al., Phys. Rev. D65 (2002) 052005. 
DO Collaboration: B. Abbott et al, Phys. Lett. B487 (2000) 264. 

[55] M. Cacciari and P. Nason, Phys. Rev. Lett. 89 (2002) 122003. 



22 




Figure 2: The differential cross section xpda /dxp for charm production in pp colhsions, (1), at 
three different laboratory energies E, for three different ways of extrapolating the gluon to very 
small X. The models are (i) a double-leading-log DGLAP extrapolation for x < 10"^ (MRST), 

(ii) a unified DGLAP/BFKL approach with an x"^ extrapolation for x < 10~^ (KMS), and 

(iii) an extrapolation with saturation effects included (GBW). Plot (d) compares the GBW 
prediction at the three energies. 
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Figure 3: The energy dependence of the relevant moment', (11), of charm production in 
proton-nucleon colhsions, as a function of the energy of the produced charm quark. The models 
are as in Fig. 2. The reason why the KMS prediction falls below the other predictions at the 
smaller values of Ec is explained in Section 3. The 7 = 2.02 moments are shown for illustration; 
in the calculation of the neutrino fluxes, the differential cross section xpda/dxp is convoluted 
with the observed primary cosmic ray flux. 
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Figure 4: The energy dependence of the relevant 'Z moment', oZ^. = J x'^'^'^{da'^/dx)dx, of 
charm production in p-air colhsions, as a function of the energy of the produced charm quark. 
The dashed curve corresponds to the GBW model (extended to include rescattering of the cc 
pair within the air nucleus). The upper dotted (B), lower dotted (gsp) and continuous {g2,p + B) 
curves respectively include the growth of the proton radius, triple-Pomeron effects and the 
combination of these two effects. The dot-dashed curve is the scaling prediction of (10), but 
with (Tinei corresponding to p-air collisions. The 7 = 2.02 moments are shown for illustration; 
in the calculation of the neutrino fluxes, the differential cross section xpda/dxp is convoluted 
with the observed primary cosmic ray flux. 
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Figure 5: The flux of prompt muon neutrinos at ground level, weigfited by E"^, for different 
clioices of tlie c cliarmed liadron fragmentation function, drf^^ /dx. Tlie curves correspond 
in descending order to assuming (i) no fragmentation dn/dx oc 5(1 — x), (ii,iii) 5{xd — 0.75xc) 
and S{xa^ — |(a^c + 1)) for Xc > 0.1, or S{xa^ — lATxc) and (iv) a Peterson et al. fragmentation 
function [43] with Ec = 0.043 [13]. In each case charm production is calculated using the GBW 
solid curve of Fig. 4. 
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Figure 6: The prompt z/^ + z/^ and v^- + fluxes, arising from cc production, fragmentation and 
decay, obtained using the three different extrapolations of the gluon to very small x (described 
in the caption to Fig. 2). For the MRST and KMS models, the A dependence is taken to be 
da/dx{p + air — i> c + . . .) = Ada / dx{j)N — c + . . .). 
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Figure 7: The prompt (a) z/^ + and (b) fr + Pj. fluxes calculated using the charm production 
cross sections corresponding to the shaded band in Fig. 4. Also shown are the conventional 
muon and tau neutrino atmospheric fluxes, where the latter originates, via neutrino mixing 
transitions, from the former. There is also a contribution to the prompt + Vr flux from 
beauty production, which is not included here, but is shown in Fig. 10(b). The prompt Ve + 
flux is equal to that for + v^, but the atmospheric flux is a factor of 10 or so less, see the 
discussion in Section 8.1. 



28 




Figure 8: The dashed and dotted curves correspond, respectively, to the prompt 



flux 



obtained by switching off the charmed hadron-air interactions and assuming that the incoming 
cosmic rays have {A) = 7 rather than A = 1. The default continuous curve (with air interactions 
and y4 = 1) is based on the GBW gluon extrapolation. In this work we assume that the c ^ Ac 



hadronization is given by xa 



for Xc > 0.1. 
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Figure 9: The differential cross section xpda/dxF for charm and beauty production in p-air 
colhsions for four different laboratory energies E. For hb production we take m?, = 4.5 GeV. 
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Figure 10: The prompt Ur + fluxes originating from cc and hh production and decay, which 
respectively arise from the Dg — > rz/^ decay, and from the S^, 5°, Bg and A;, semileptonic rv^ 
decay modes. The upper plot shows the breakdown into the direct v-r contribution (continuous 
curves) and the indirect r ^ v-r contribution (dashed curves). The lower plot shows the total 
prompt Vj. + Vj. flux, together with its components of cc and hb origin. Also shown is the non- 
prompt Vt- + Vj. flux arising from oscillations from the conventional atmospheric v^J^ 
flux. 
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Figure 11: The continuous curve is the prompt z/^ + z/^ flux from horizontally- incident cosmic 
rays. The dashed curve, which corresponds to vertically-incident cosmic rays, is taken from 
Fig. 7. 
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